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ABSTRACT 

Numerical simulations of quantum chromodynamics at nonzero temperature pro- 
vide information from first principles about the physical properties of the quark 
gluon plasma. Because the lattice approximation can be refined indefinitely, 
results of lattice simulations now provide the most reliable basis for our under- 
standing of the nonperturbative characteristics of the plasma and of the high 
temperature phase transition. Following a brief overview of the methodology of 
lattice gauge theory at nonzero temperature, recent results and insights from 
lattice simulations are discussed. These include our understanding of the phase 
diagram of QCD, the nature of the phase transition, and the structure of the 
plasma. 
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1. Introduction 

In 1974 Wilson proposed a formulation of quantum chromodynamics (QCD) on a 
discrete space-time lattice.^ A few years later Creutz demonstrated that this formu- 
lation was a promising basis for the successful numerical simulation of the theory.^ 
The simulation method is suitable for studies of QCD at nonzero temperature in 
thermal equilibrium and, at least at present, at zero baryon chemical potential. It is 
the most promising currently available method for deducing nonperturbative charac- 
teristics of QCD at zero or nonzero temperature, starting from first principles. Given 
suitable algorithms and computing power, the only approximations, namely, finite 
lattice spacing, finite volume, and unphysically large quark masses, can be refined 
indefinitely. Thus lattice simulations provide the basis for much of our theoretical 
understanding of the characteristics of QCD at nonzero temperature. 

Are the predictions of lattice QCD relevant to the experimental effort? The vi- 
olently expanding product of a high energy heavy ion collision can scarcely be ap- 
proximated as an equilibrium plasma at zero baryon chemical potential. Indeed, even 
a hydrodynamical description that assumes local equilibrium could be questionable. 
Thus despite our confidence in lattice QCD, it is still essential that we formulate good 
models to bridge the gap between the relatively solid predictions of the computer sim- 
ulations and the experimental environment. 

In the past few years lattice simulations have achieved a reasonably consistent 
qualitative picture of the high temperature behavior of QCD. What has emerged is 
a far more subtle and intriguing picture of the quark gluon plasma and the nature 
of the phase transition than earlier caricatures would have suggested. At very high 



temperatures we have a plasma that behaves in bulk roughly like a free gas of quarks 
and gluons, yet it retains features of confinement that are evident in long-range cor- 
relations. Although there may be no phase transition, there is at least a dramatic 
crossover accompanied by dramatic changes in the structure and symmetries of the 
plasma. In this chapter we will be taking stock of what lattice QCD tells us about 
hadronic matter in thermal equilibrium. Particular emphasis will be placed on de- 
velopments since the last volume, which featured an excellent review by Karsch."^ For 
the past several years, the annual "Lattice" conference has regularly featured sessions 
on nonzero temperature QCD. Recent proceedings can be found in Refs. 4,5,6,7. 

2. Nonzero Temperature QCD on the Lattice 

Here we give a very brief sketch of the basis for lattice gauge theory at zero and 
nonzero temperature. Our purpose is to indicate how the Euclidean space-time lattice 
with periodic boundary conditions arises naturally in the quantum Gibbs ensemble, 
and to show how the temperature of the ensemble determines the imaginary time 
dimension of the lattice. For a more complete exposition, there are some excellent 
texts.® 

2.1. Feynman Path Integral from the Gibbs Ensemble: A Thumbnail Sketch 

Nonzero temperature simulations of a field theory are designed to calculate oper- 
ator expectation values on the quantum Gibbs ensemble at temperature T: 

{O) = Tr((9e-^/^)/Tre-^/^, (1) 

where H is the hamiltonian for the field theory, O is any operator, the trace is over 
all physical states in the Hilbert space, and we use units in which the Boltzmann 
constant k is one. The operator exp[—H/T) is just the standard time evolution 
operator exp[—iHt) evaluated at an imaginary time 

t = -i/T. (2) 

This expression is the origin of the imaginary time coordinate and the important rela- 
tionship between the temperature and the extent of the Fuclidean spacetime volume 
in imaginary time. 

Several steps are required in order to convert the trace over states in the Gibbs 
ensemble (1) into a multidimensional integral over lattice variables. Let us consider 
briefiy how this is done for a scalar field theory based on a single field (p. First the 
hamiltonian H is formulated on a three-dimensional lattice with lattice constant a. 
The continuous field then takes on values ^(x) on each of the sites x of a three- 
dimensional lattice. The trace over states is written on a complete orthonormal basis 



|^(x)) in which the field is diagonaL To facihtate the estimation of the time evolution 
operator, the imaginary time interval [0, l/T] is then subdivided into Nf steps for large 
Nf^ and the trace (partition function) is rewritten in the form 

Z = Tre-^/^ = Tr(e-^/^'^e-^/^'^ . . . e"^/^'^). (3) 

This step produces the lattice discretization in imaginary time. Each of the infinitesi- 
mal time evolution operators (transfer matrices) in the product is written as a matrix 
on the field-diagonal basis. This is done through the completeness relation 

1= /n#(x)|<^(x))(<^(x)|. (4) 

Since this relation is inserted between each factor, it is convenient to introduce an 
extra label r to distinguish the multiplte integration variables: ^(x, r). The extra 
variable is naturally taken to be a discrete imaginary time variable, leading to a 
classical field variable ^(x, r) defined on a four-dimensional Euclidean lattice. The 
discrete time values are 

T = a^k for A; = 0,1,2,. .. ,iVt - 1. (5) 

where = 1 /N^T is taken to be the lattice constant in the imaginary time direction. 
In terms of this labeling the partition function then becomes 

Z = /n#(x,r)(<^(x,0)|e-^/^'^|<^(x,(iV,-l)aO) (6) 

X,T 

((^(x, (N, - l)aO|e-^/^'^|<^(x, (N, - 2)a,)) . . . 

(<^(x, 2a0|e-^/^'^|<^(x, aO) (<^(x, aO|e-^/^'^|<^(x, 0)). 

Since we are taking a trace, we have built in the requirement that ^(x, Nfat) = ^(x, 0). 
This is the origin of periodicity in imaginary time. An explicit evaluation of the 
transfer matrix elements leads to Eeynman's remarkable path integral formula for the 
quantum partition function 

Z(r) = Tre-^/^ = fl[d<P{x)exp[-S{<P,T)], (7) 

where S{(f>^ T) is the imaginary time classical action for the field configuration cf>{x) for 
X = (x,t) on the Euclidean space-time lattice of dimension X Nf. The integration 
is over all possible choices of the field values on the lattice. 

Any observable (9 is a function of the field cf). The expectation value (1) is then 



{O) = / Y{d<i>{x)0{<i>)eM-S{<i>.T)]lZ{T) 



(8) 



2.2. Monte Carlo Methods 



As sketched in the previous section, the Feynman path integral formulation re- 
duces quantum statistical mechanics to an integration over classical variables. For- 
mulated on a lattice, the problem is reduced to a multidimensional integration. For a 
wide class of actions, the weight factor exp(— S*) of the integration is positive definite, 
and the integration can be done effectively using a variety of Monte Carlo sam- 
pling methods, using the weight factor as a probability. Among these are heat bath. 
Metropolis, and molecular dynamics methods. The basic idea is to produce a large 
biased sample of points {(f>i[x)^ z = f , . . . , A^conf } in the space of the multidimensional 
integration. These are commonly called "configurations", and are characterised by 
choosing a particular value of the field on each space-time point. The sample is biased 
so that the probability of encountering a configuration is proportional to the weight 
factor 

Pmx)})cxexpi-S). (9) 

On such a biased sample, the operator expectation value is simply the average of its 
values on the sample configurations: 

((9)= hm ^E'wO- (10) 

The algorithmic challenge lies in formulating an efficient sampling method that pro- 
duces a desired variance at the lowest computational cost. 



2.3. Lattice QCD 



2.3.1. Pure Yang-Mills Theory 

The Feynman path integral is expressed in terms of the action in a Fuclidean space 
time. The Wilson formulation of the Fuclidean action for quantum chromo dynamics 
starts from a regular hypercubic lattice with equal space and time lattice constants 
a = a^. The gauge vector potential A°'^{x) defines the gauge connection 

[/,,^ = exp(z(/aA^(x)A,/2) (ff) 

between the site x and the nearest neighbor site x -\- jj.. (The lattice vector of length 
a in the // direction is /i.) Here g is the gauge coupling constant and Ac are the 
usual generators of the SU(3) Lie algebra. This SU(3) matrix is called the gauge link 
matrix. There is one such variable for each of the links connecting nearest neighbors 
on the lattice. A forward connection for a given link is associated with the matrix U 



and a backward connection for the same link is associated with the adjoint of that 
matrix 

The plaquette variable is defined on a unit square on the lattice as the product of 
the connections around the square. 

□ = TtU,,{x) = TtU,{x)U,{x + jl)Ul{x + p)Ul{x) (12) 

The trace of any such product of gauge connections around a closed path is gauge 
invariant. The plaquette variable is related to the SU(3) color Maxwell field strength 
in the continuum limit. 

limReTr[/,4x)/3 = 1 - ^[F^,,{x)]' (13) 
The continuum Euclidean action for a pure gluon field is 

S, = J d'xl[F^Jx)]\ 
Wilson suggested the lattice approximation 

= E E 6//[ReTr[/,.(x)/3 - 1] (15) 

Creutz first demonstrated the feasibility of using this simple approximation in 
numerical simulations. It has served as the basis for a great many studies since then. 
Improvements, some very promising, have been proposed. They add terms to the 
action formed from products of the gauge connection around larger loops.^'^°'^^ The 
objective is to remove lattice artifacts as nearly as possible and facilitate the approach 
to the continuum limit. 

2. 3. 2. Including Quarks 

Incorporating fermions into the functional integral demands an additional effort. 
The Pauli exclusion principle requires that they be introduced as anticommuting 
Grassmann numbers rather than the ordinary commuting numbers of the boson fields. 
Since we compute with ordinary numbers, it is then necessary to complete the inte- 
gration over the fermion degrees of freedom by hand. Fortunately, this is easy to do. 
Unfortunately, the resulting nonlocal effective action vastly increases the computa- 
tional cost. The result, however, is a simulation of full QCD. 

An additional complication with fermions is a difficulty in controlling the number 
of fermion species. We would like to construct a theory that reproduces faithfully 
the chiral symmetry of the continuum theory at zero quark mass. The lattice regu- 
larization forces us to make a difficult choice. Either we give up chiral symmetry or 
we must have a doubling (usually a few redoublings) of the number of quark species. 



(14) 



There are two popular lattice fermion formulations corresponding to these choices. 
One is called the staggered or Kogut-Susskind fermion formulation^^ and the other 
the Wilson fermion formulation.^ The hope is that in the continuum limit, either 
choice takes us to the one and only continuum theory. 

For present purposes we merely write down the lattice fermion actions for these two 
choices. For further details the reader should consult Ref. 8. In the Wilson fermion 
formulation the quark field for each fiavor is represented as a Dirac color spinor V'j(^) 
on each lattice site x with a three- component color index c and a four- component 
Dirac spin index j . The fermion action is then 

'^7 = J2{ - Kj2[i^ix + + 7A*)V'(a^) + i^ix)ir - 7^)V'(a^ + A)] > (16) 

where k, = l/(2am -|- 8r) and r is usually taken to be 1. The summation over color 
and spin degrees of freedom is implicit. With r = the fermion action is chirally 
symmetric at zero quark mass and describes 16 degenerate fermion species. With 
r ^ the degeneracy is lifted at the expense of destroying chiral symmetry in the 
zero quark mass limit. 

In the "staggered" fermion formulation the quark field is represented as a color 
spinor t/j'^[x) with no explicit Dirac spin degree of freedom. In effect, four spin and 
four fiavor components are distributed over each hypercube of dimension 2"*. The 
fermion action is 

S'j = ^ < 2am%l){x)%l){x) ^ a^{x)['il){x)%l){x + jl) - 4){x jj)%l){x)] \ (17) 

The phase factors a^(a;) are diagonalized Dirac matrices. The summation over the 
color index is implicit. The theory is chirally symmetric at zero quark mass, but 
there are four degenerate fiavors. Such a fiavor symmetry is unnatural, but there are 
methods for reducing the effective fiavor number.® Because fiavor rotations involve 
fields on different lattice sites, except for U(l) transformations, which are diagonal 
in space, fiavor rotations are restricted to a discrete lattice subgroup. Thus instead 
of the full SU(4) xSU(4) chiral symmetry, the symmetry consists of a U(l)xU(l) 
subgroup plus a discrete subgroup. Since a restoration of rotational invariance in the 
continuum limit is anticipated, it is also expected that the full chiral symmetry will 
be restored simultaneously. On a coarse lattice members of the same fiavor multiplet 
are not necessarily degenerate. For example, it is popular to measure the masses of 
two members of the pion multiplet, often called tt and 7r2, with local operators in the 
staggered fermion scheme. The mass ratio m^/m^2 serves as an indicator of progress 
toward the continuum limit, since it should approach one. 

The total action in either case is the sum of the gauge and fermion parts 



S — S g ~\~ S j . 



(18) 



The partition function is then given by 



Z = j dUdipdipe-'^ (19) 

As we have noted, the fermion fields must be integrated out explicitly. In either 
formulation, the fermion fields enter in a bilinear form, 

Sj = y^M{U)tP (20) 

where the fermion matrix M[U) is a matrix with row and column indices labeled by 
the spatial coordinate as well as color and, if necessary, spin. Integrating out the 
quark degrees of freedom leads to 

/ n dtlj{x)diP{x) exp[tljM{U)ip] = det M{U) (21) 

Thus the effective gauge action is 

Ses{U) = S,{U) + log det M{U) (22) 
and the partition function is then 

Z = [ J[/e-^=ff(^) (23) 



The second term in the effective action is the fermion determinant. It depends on 
the gauge fields and induces a nonlocal gauge field interaction. This feature vastly 
increases the computational effort. It is beyond the scope of this review to dis- 
cuss the various methods for accommodating the fermion determinant in a numerical 
simulation.® The most effective method uses a molecular dynamics approach. With 
tricks it is possible to approximate the induced effect upon the gauge field of any 
number of fermion species. The more elegant "exact" simulations require four species 
of staggered fermions or two species of Wilson fermions. 

2.3.3. Asymptotic Freedom 

Since our goal in lattice QCD is to regulate and approximate the continuum 
theory, it is crucial that there be a meaningful lattice continuum limit. Fortunately, 
QCD (and lattice QCD) is an asymptotically free field theory with an ultraviolet 
fixed point. As the lattice coupling g is decreased, the lattice constant a, measured 
in physical units, also decreases. For small enough </, it decreases according to the 
perturbative scaling relation 

aA = e-i/(2/^«^')(/3o5')-''^/('''o) (24) 



where /3o = (ll-2iVj/3)/(167r2) and /3i = (102-38iVj/3)/(167r2)2 and A sets the scale. 
It is determined in lattice regularization from a physical quantity, such as the proton 
mass. A proton may occupy only a few lattice sites at large </, but as g is decreased, 
the proton occupies more lattice sites, so is represented with increasing resolution. 
It follows from the relationship between g and lattice scale a, that increasing 6/g^ 
corresponds to increasing the temperature T = l/[Nta). Thus by varying 6/(/^, it 
is possible to study a range of temperatures without changing the number of lattice 
points. 

2. 3. 4. Quenched Approximation 

The "quenched" approximation to QCD amounts to carrying out the simulation 
with the fermion determinant set to one, cutting the computational effort by orders 
of magnitude. In this case the thermal ensemble consists only of gluons. Quite useful 
results can be obtained in this way at zero temperature, where it has been difficult to 
find cases where vacuum fermion loops seem to make a difference. At nonzero tem- 
perature, however, fermions do make a significant difference in the thermal ensemble. 
The phase structure changes when quarks are included. Thus the quenched approxi- 
mation at nonzero temperature may give suggestive results, but fermions cannot be 
neglected if precise contact with nature is needed. 

2. 3. 5. Nonzero Chemical Potential 

The partition function for full QCD, Eq. (19), represents the grand canonical 
ensemble at zero baryon chemical potential. The fermion determinant is real and 
nonnegative for two fiavors of Wilson fermions or four fiavors of staggered fermions. 
Thus the factor exp(— S'efr) is a suitable probability weight for a Monte Carlo sim- 
ulation. It would be very useful to be able to carry out simulations at nonzero 
chemical potential. The problem is obviously important, since the debris of a heavy 
ion collision necessarily includes regions of nonzero chemical potential and the cores of 
neutron stars are obviously baryon rich. Unfortunately at nonzero chemical potential, 
we encounter a fundamental technical problem: The fermion determinant becomes 
complex, spoiling the probability weight. Various tricks have been proposed to evade 
this problem. They include carrying out simulations in the quenched approximation 
or at zero chemical potential, and incorporating a correction factor to simulate the 
effects of a nonzero chemical potential. Such methods are limited to a very small 
chemical potential or a very small volume. Thus they may succeed in giving hints 
about what happens at a small chemical potential, but none has been successful in 
achieving a result that can be confidently extended to the thermodynamic limit of 
inifinite volume. Thus this problem remains an open challenge. For a recent review. 



see Ref. 13. 



2.4. Lattice Ohservahles 

2.4.1. Quark Propagator 

To determine properties of the thermal ensemble one measures a variety of ob- 
servables. Many of the important observables involve quarks. Thus, for example, we 
need the quark propagator 

Wz)y^{y)) = I Y{dU,{x)d4^{x)d4^{x)e-'4^{z)4^{y) (25) 

The integration over quark variables is again carried out explicitly, giving 

m^YM) = j Y{dn,{x)M-\z,y,U)eyi^{-S,s{U)) (26) 

or the gauge field average of the inverse of the fermion matrix. (We have suppressed 
the color and spin labels for simplicity.) In a more compact notation, we may write 

{i^{z)i,{y)) = {M-\z,y,U))u (27) 

where the subscript U in the expectation value on the right side indicates that it is 
taken with respect to the effective gauge action. 

2.4.2. Hadron Propagator 

If an observable is not gauge invariant, the integration over gauge variables gives 
zero. The quark propagator we have been discussing can be defined in a specific 
gauge. A hadron propagator, on the other hand, is gauge invariant. For example 
the operator il){x)Vil){x) is an interpolating field for a quark-antiquark meson with a 
particular Dirac matrix V determining the spin and parity. Thus a meson propagator 
can be extracted from the correlation 

{i^{z)vi,{z)i,{y)vi,{y))u - {mm^))u{mm^))u 

= {TT[M-\z,y,U)VM-\y,z,U)V])u 

- {TT[M-\z,z,U)V]TT[M-\y,y,U)V])u (28) 

- (Tr[M-i(0, 0, [/)r])c/(Tr[M-i(0, 0, U)V])u 

The second term on the right side makes a contribution only for fiavor singlet mesons. 
It represents a coupling to gluon intermediate states. The last term, the vacuum 



disconnected term, on either side contributes only for flavor singlet mesons with 
vacuum quantum numbers, such as the chiral condensate order parameter t/jt/j. 

In this way all hadron propagators are constructed from products of the inverse of 
the Dirac matrix, averaged over gauge held conflgurations with a weight determined 
by the effective gauge action. 

2.4.3.Polyakov Loop 

An important observable simulates the effect of introducing a static spinless quark 
into the ensemble. It is formed from the product of a string of gauge link matrices 
along a line in the time direction. The observable that introduces a static quark at 
site X is 

Nt-l 

P(x) = Tr n t^4(x,T) (29) 

T = 

where the trace is over the color degrees of freedom. The static quark world line is 
closed by virtue of the periodicity of the lattice. A static antiquark is introduced by 
the complex conjugate variable. 

The change in the free energy of the ensemble caused by the addition of a single 
heavy quark is given by 

exp[-/(r,m,)/r] = (P(0)/3)c/ (30) 



2. 4. 4. Heavy Quark Potential 
The thermal heavy quark potential is determined from 

exp[-y(r, T)/T] = {P{0)P\r)/9)u (31) 

The quantity V{r^ T) is more precisely the change in the free energy of the ensemble 
caused by adding a spinless quark and antiquark pair at separation r. 

3. Phase Structure 

In nature the quark masses assume their physical values, of course. In numerical 
simulations, however, it is possible to adjust quark masses and other parameters 
to gain more insight into the phase structure of QCD. Thus a phase diagram can be 
constructed in the multidimensional space of the temperature T, the quark masses nii 
for z = 1, . . . ,Nf flavors and the corresponding chemical potentials IJ4. The majority 
of simulations have been carried out at zero chemical potential with two or four 
flavors of equal mass quarks. However, there are a few simulations that approach 
a more physical quark mass spectrum with two equal mass light quark flavors and 
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Figure 1: Schematic phase diagram in quark mass and temperature for QCD. 

one strange quark. Because simulations become more costly as the quark mass is 
decreased, simulations are done at an unphysically large quark mass, requiring an 
extrapolation to smaller physical values. 

At infinite mass, quarks are absent from the thermal ensemble, and QCD be- 
comes a pure SU(3) Yang-Mills theory. It is well established that this theory has a 
first order nonzero temperature phase transition. Indeed the attendant coexistence 
of plasma and confined phases at the critical temperature inspired speculation about 
phase boundaries and bubbles in plasma production and cooling. At the opposite 
extreme of zero quark mass, QCD becomes a chirally symmetric theory. The sym- 
metry is spontaneously broken at low temperature, leading to a zero mass pi meson. 
As temperature is increased, a phase transition occurs, leading to restoration of the 
broken symmetry and the appearance of chiral multiplets in the hadron spectrum. 
Thus we have the schematic phase diagram for QCD as a function of a single SU(Aj) 
fiavor symmetric quark mass and temperature shown in Fig. 1. The first order con- 
finement /deconfinement phase boundary has been extended to finite quark mass. 
Whether the chiral phase boundary extends to nonzero quark mass as shown or oc- 
curs only at zero mass depends on the number of fiavors. With three or four light 
quarks quite general arguments can be made in favor of a first order chiral phase 
transition.^'* In that case an extended chiral phase boundary is expected. But with 
two light fiavors a second order phase transition is not ruled out,*"* in which case it 
may occur only at zero quark mass. Thus symmetry and universality considerations 
do not require a phase transition at physical values of the quark masses. To establish 
its existence, we turn to numerical simulation. 
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Figure 2: QCD phase diagram for 2+1 flavors as a function of a degenerate up and 
down quark mass niu^d = = and a strange quark mass m^. The heavy hne is 
the critical phase boundary. 

Are these expectations conflrmed in numerical simulations? For four flavors of 
equal mass quarks numerical evidence suggests that the chiral phase boundary extends 
to nonzero quark mass, but it is likely that a gap separates the chiral and deconflning 
phase boundaries as sketched in Fig. 1."^'^^ For two flavors of equal mass quarks there 
appears to be no phase transition at a small, but nonzero quark mass.^^'^"* For two 
equal light quark masses and one strange quark, we view the phase diagram from a 
different perspective in Fig. 2.^'' Such a phase diagram is motivated by an analysis 
of the corresponding sigma models in mean held theory, augmented by an analysis of 
quantum fluctuations. ^'^'^^'^^ The diagram indicates for which range of the two mass 
parameters a nonzero temperature chiral phase transition occurs and whether the 
phase transition is expected to be flrst or second order. When the strange quark 
mass is large, the thermal ensemble is effectively a two-flavor ensemble, and we have 
a second order phase transition only at zero up and down quark mass. However, 
when the strange quark is sufficiently light, we recover the first order chiral phase 
transition expected of the three-fiavor ensemble, which extends to nonzero quark 
mass as sketched. Whether QCD conforms to this expectation in detail remains to 
be established. 

So is there a phase transition at physical values of the quark masses? Simula- 
tions with 2 + 1 fiavors in the staggered and Wilson fermion schemes both support 



the existence of the first order region sketched in Fig. 2, but do not agree on its 
extent. Staggered fermion simulations of the Columbia group found a first or- 
der signal for {mu^d^'m-s) ~ (15,15) MeV, but none for {niu^^rris) (15,30) MeV, 
suggesting no phase transition at physical values. A recent Wilson fermion simula- 
tion by the Tsukuba group found a first order signal with quark masses as large as 
{mu,d,ms) ^ (140,140) MeV and {mu,d,ms) ^ (0,400) MeV,^°'^^ allowing a phase 
transition at physical values. Although the staggered fermion approach for studying 
nonzero temperature QCD is more mature than the Wilson approach, until algorith- 
mic improvements lead to consistency, we cannot be certain whether there is a phase 
transition. Nonetheless, it is likely that if there is a phase transition at all, it is weak. 
Still, as we shall see, even if there is no phase transition, there is a dramatic crossover. 
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Figure 3: Polyakov loop susceptibility vs. from Karsch and Laermann. Solid 

symbols are directly simulated. Open symbols are derived by reweighting. Peaks 
locate the crossover. 

4. Temperature of the Phase Transition 

As we have seen, light quarks make a significant difference in the character of the 
phase transition, so cannot be neglected in a realistic study of thermodynamics. The 
temperature of the phase transition or crossover in the presence of dynamical quarks 
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Figure 4: Polyakov loop vs. temperature from Ref. 27. The line segments in- 
dicate the range of uncertainty in locating the maximum slope. Data are from 
Refs. 24,26,30,37,26. 

is also an important indicator of convergence to the continuum limit. Knowing its 
approximate value as well is of considerable phenomenological importance. 

There are two methods to locate the crossover. The most elegant locates the peak 
in a susceptibility, such as that based on the Polyakov loop."^° 

XL = NM^ePY) - {RePn (32) 

Combined with a Ferrenberg-Swendsen reweighting analysis to interpolate between 
simulation points, the method permits a clean determination of the crossover coupling, 
as shown in Fig. 3. A second method locates the peak in the derivative of a sensitive 
observable, as illustrated for the Polyakov loop variable in Fig. 4. 

The equilibrium temperature in a lattice simulation is the reciprocal of the lat- 
tice extent in the imaginary time dimension, or l/Nfa. To express the crossover 
temperature in physical units requires the measurement of an additional quantity at 
the same lattice parameters where the crossover is observed, for example, the (zero 
temperature) rho meson mass, also determined in dimensionless lattice units as nipa. 



0.3 



0,2 




0.1 



Wilson 
Staggered 
Wilson 
Staggered 
Wilson 
Staggered 
* N^=12 Staggered 



o Nt=4 
+ Nt=4 
□ Nt=6 
K Nt=6 
Nt=8 
X Nt=8 



0.0 



0.0 



0.2 0.4 0.6 0.8 

(111^2/111^)^ 



1.0 



Figure 5: From Ref. 24. Temperature of the crossover in units of the rho meson mass 
with two hght quarks vs. the squared ratio of the pion to rho mass. For points from 
the staggered simulations, the mass of the second (non-Goldstone) pion is used. The 
curved hne segment shows the error bar for the new = 8 staggered point. The 
vertical dashed line indicates the physical mass ratio. Data are from Refs. 24,26,37,50 

Thus 

Tjnip = l/{Ntamp) (33) 

To set the scale in this way requires a separate simulation at low temperature. From 
the experimental value the temperature is then converted to physical units: 

r,(MeVp) = 770r,/m, = 770/(iVtam,) (34) 

The result must be interpreted cautiously, however, since in present simulations with 
dynamical quarks, the ratio of the nucleon to rho meson mass is unphysically high. 
Setting the scale with the nucleon mass instead would therefore result in a still lower 
temperature. There is some hope for future improvement: Through a series of ex- 
trapolations it has been possible to reach physical values of the nucleon to rho mass 
ratio in recent zero temperature quenched calculations for hadronic states constructed 
from both Wilson and staggered valence fermions.^"^ The same accomplishment with 
dynamical quarks included will take considerably more effort. Another problem af- 
fecting the determination of the crossover temperature is the need for a rather long 



extrapolation to physical quark masses. In present simulations with dynamical quarks 
the zero temperature mass ratio m^/mp lies in the range 1/3 to 1/2, twice to three 
times the experimental value of 0.18. This ratio can be reduced to its physical value 
by reducing the quark mass considerably. So for now we use the notation MeVp to 
call attention to the assumptions made.^^ 

The temperatures thus determined are plotted as a function of the square of the 
ratio of the pi to rho mass in Fig. 5.^"* This figure includes results from simulations 
in the Wilson as well as staggered fermion schemes. As we have noted, the Wilson 
scheme breaks chiral symmetry explicitly and the staggered fermion scheme breaks 
the continuous fiavor symmetry explicitly. The breaking of these symmetries results 
in an unphysical value for the chosen ratio. Both symmetries are expected to be 
restored in the continuum limit, permitting approach to the physical value. Thus 
progress toward the continuum limit is measured by movement toward the physical 
ratio. There is an encouraging consistency in the results, which place the crossover 
temperature at approximately 

r, = 140 - 160 MeVp. (35) 

By contrast, in a gluon ensemble without quarks (quenched approximation), the tem- 
perature of the phase transition is about 260 MeV.^^ For further details concerning 
the construction of the temperature scale, please see Refs. 26,27. 

5. Consequences of Chiral Symmetry Restoration 

The QCD phase transition or crossover inherits characteristics of both the decon- 
finement and chiral phase transitions. An immediate consequence of chiral symmetry 
restoration is the vanishing of the chiral order parameter (V'V') iii the limit of zero 
quark mass. This effect is observed. Moreover, the vanishing is quite abrupt. We 
begin this section with a discussion of these results. The appearance of critical behav- 
ior is a second consequence of chiral symmetry restoration, depending sensitively on 
the values of the up, down, and strange quark masses. Critical behavior is expected 
along the curve shown in Fig. 2. However, since so little work has yet been done 
with two plus one fiavors, we restrict our discussion of critical behavior to the case of 
zero quark mass in two-fiavor simulations. Finally, the restoration of chiral symmetry 
also produces dramatic shifts in the hadron spectrum, leading to the emergence of 
chiral multiplets. These multiplets appear in long range screening correlators, e.g. 
the high temperature Yukawa interaction between two nucleons and are readily ob- 
served in numerical simulations. They do not appear to have any directly observable 



consequences, but their existence imposes constraints on phenomenological models. 



5.1. Chiral Order Parameter 

The dramatic fallofF in the chiral order parameter (V'V') has long been used to 
locate the high temperature crossover. Figure 6 from Boyd et al.^^ gives a compilation 
of results for the ratio {t/jt/j)[T) / {t/jt/j)[T = 0). This figure collects results for a variety 
of staggered quark fiavors, including quenched [Nj = 0) at Nf = 4^^ and Nf = 8,^® 
two fiavors 



anir 



0.02, three fiavors (two light, one heavier) with amu,d = 0.0125 
and anig = 0.25,"^^ and four fiavors. "^^ The quenched results are extrapolated to zero 
quark mass. The others are at small quark mass, but not extrapolated to zero. The 
zero temperature value is from Ref. 33. 

Also evident in Fig. 6 is a remarkable constancy in the order parameter for T < 
0.92Tc. The same study also found no significant variation in the pion decay constant 
over this temperature range. ^® Thus any experimentally observed shifts in hadron 
properties would indicate temperatures very close to or above the critical temperature. 
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Figure 6: From Ref. 28. Chiral condensate normalized to its zero temperature value 
as a function of T/Tc. Zero fiavor results are extrapolated to zero quark mass. The 
others are not. See text for details. 



5.2. Critical Behavior 

If QCD falls into the proposed universality class of sigma models with a second 
order chiral phase transition at the boundaries indicated in Fig. 2, certain scaling 
relations apply near the critical point. "^"^ The question is crucial for a successful ex- 
trapolation to small quark mass and large volume near the phase transition. A correct 
extrapolation requires the correct critical exponents. Thus an important test of the 
proposed phase structure is to determine whether QCD has the expected critical 
behavior. To emphasize that the conclusion is not foregone, a recent analysis of the 
three-dimensional Gross-Neveu model by Kocic and Kogut questions the conventional 
wisdom that places QCD with its composite scalar mesons in the same universality 
class as sigma models with their elementary scalar mesons. "^^ Instead, in a detailed 
simulation Kocic and Kogut found mean-field scaling. 
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Figure 7: Crossover coupling v.s. quark mass for two fiavors of staggered quarks from 
Ref. 30. Solid line is a fit with 0(4) critical exponents. Dashed line is a free fit. 

Let us focus on the two fiavor theory in the staggered fermion scheme, corre- 
sponding to the upper portion of the nig axis in Fig. 2. In this case 0(4) universality 
is expected in the continuum limit. To be more precise, in the staggered fermion 
scheme, one expects 0(2) critical behavior on coarse lattices where the fiavor sym- 



metry breaking of the staggered scheme is significant, and 0(4) as the continuum 
hmit is approached. Recent work by Karsch and Karsch and Laermann attempts a 
determination of some of the critical exponents of QCD."^°'"^^ They exploit the stan- 
dard sigma model analogy between QCD and a magnetic system. In this analogy the 
quark mass plays the role of an external magnetic field and {^^tp^ plays the role of 
the magnetization. 

The standard analysis of critical behavior begins with the assumed scaling of the 
critical contribution to the free energy in a magnetic system. "^"^ This contribution is 
singular and dominant at small field close to the zero field critical temperature Tc{Q). 
The scaling property, expressed in terms of the scaled temperature t. 



t = [T- r,(o)]/r,(o) (36) 



and magnetic field h is 



Ut{t,h) = h-^U,{hyH,hy'^h). (37) 

From this scaling behavior one derives a scaling relation for the critical contribution 
to the magnetization s = —dfcrit/dh: 

s{t,h) = h^l^y{x) (38) 

where 

X = th^l^^ (39) 

and y{x) is a scaling function. Here fi and 8 are the critical exponents appropriate 
to the universality class. An important consequence of this result is that a crossover 
peak in the susceptibility Xh = ds/dh moves along a curve of constant x = Xpc as h 
and t are varied. Thus if critical scaling holds, once the pseudocritical temperature 
is known at one /i, it can be predicted at all h. 

In QCD the quark mass plays the role of the magnetic field and (^tptl)"^, the mag- 
netization. Specifically, Karsch suggests the identification 

h = mjT = anigNt (40) 
t = -6/gliO,N,), (41) 

where 5'c(0, N^) is the critical gauge coupling at zero quark mass for a particular N^.^^ 
Scaling then predicts how the critical gauge coupling chang function of the 

quark mass: 

Q/glM = 6/gliO) + imjrf^'. (42) 

Using this expression Karsch presented an analysis of the crossover for Nf = 4, 6, 
and 8 for two quark fiavors."^^ The agreement is quite encouraging. The addition of 
new data at = A permits a more refined analysis, shown in Figure 7."^° Their best 
fit critical exponent l/jSS is 0.77 it 0.14, in slight disagreement with the 0(4) value 
0.55(2), but consistent with the 0(2) value 0.60(1) and the mean field value 0.67. 



Karsch and Laermann also introduced a new cumulant 



^ ^ d\n{i)i){%lg^,mq)) ^ 1 _ xy'{x) 
dlnniq 8 l36y[x) 

that evaluates the critical exponent A=l/6atx = t = 0. They obtain 0.21 < 1/6 < 
0.26 consistent with the 0(4) value 0.208(2) and 0(2) value 0.2080(3), and somewhat 
inconsistent with the mean field value 0.33. Thus the Kocic-Kogut scenario cannot 
be decisively excluded. 

To push these results further we can test the full scaling relation (38) in simulations 
with two fiavors of staggered fermions over the wide range of currently available 
Nf. For this purpose we use slightly different variables to permit comparison among 
different and to avoid quantities with anomalous dimensions, namely, 

h = ml{m,,T = 0)/ml{m,,T = 0) (44) 

t = [r-r,(o)]/r,(o) (45) 

s = h-^m{i}il){mq,T))lT^ (46) 

The scaling relation (38) then gives a universal function 

y{x) = h-^-^l^mq{4)il){m,T))lT* (47) 

with X = th^l^^ . The extra factor is needed to compensate for the quark mass 
factor niq. 

To implement the alternate variables (46) one must know the zero temperature pi 
and rho masses over the parameter range of the nonzero temperature analysis. This 
is done by constructing an interpolation of known hadron masses.^^'^'' 

This analysis was applied to data for {^V'V') from several groups."^'' Setting the 
critical exponents 8 and fi to their 0(4) values and adjusting the sole remaining free 
parameter Tc{Q) in 10 MeV increments to get the best agreement gives the result 
shown in Fig. 8 for T^i^) = 150 MeV. With the exception of the Nf = 12 data, 
the scaling agreement is rather good. At this level it is not possible to distinguish 
0(4) from 0(2) and mean field critical behavior. Setting Tc{0) to 140 MeV or 170 
MeV worsens the agreement noticeably, but 160 MeV gives comparable consistency. 
Obviously a host of systematic errors, including finite volume effects and deviations 
from continuum scaling enter the analysis, so refinements are certainly needed before 
the method can serve as a definitive test of critical behavior. The most glaring 
inconsistency in this figure comes from the preliminary = 12 data. Increasing 
Tc{0) for only this data set to 160 MeV and plotting it with other data computed 
for Tc{0) = 150 MeV brings the Nf = 12 data at the lower quark mass into good 
agreement. Thus the discrepancy could be caused either by a gradual upward shift 
in the crossover temperature as the lattice spacing is decreased, or by an erroneous 
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Figure 8: Scaled (V'V') ^•5- scaled temperature with 0(4) critical exponents. Data are 
from Refs. 24,26,30,37. 

extrapolation of the hadron spectrum above 6/(/^ = 5.7, or by a failure of the scaling 
hypothesis over this parameter range. 

5.3. Screening Spectrum 

When a symmetry is restored in a hamiltonian system, we expect eigenstates 
to organize themselves into multiplets belonging to representations of the symme- 
try group. The chiral group of QCD is U^(l) X SU(Nj) X SU(Nj), depending on the 
number of degenerate quark flavors N j. At low temperature the flavor singlet U^(f) 
symmetry is explicitly broken by the gauge anomaly and the flavor nonsinglet symme- 
try SU(N/)xSU(N/) is spontaneously broken to the flavor group SU(N/). The gauge 
anomaly receives contributions from instanton gauge configurations, which are ex- 
pected to be suppressed at high temperature."^® The spontaneously broken nonsinglet 
symmetry is also expected to be restored at high temperature. If chiral symmetry 
is restored, which part of the symmetry is restored? The answer to this ques- 
tion helps us understand what soft modes may be present at the crossover. Critical 
fiuctuations drive soft modes, and may give rise to collective behavior with observ- 
able consequences.'*^ We may seek the answer by examining the spectral multiplets. 
For example, under SU(2) X SU(2) the expected low-lying meson multiplets are the 



quartets {vr,(T} and {?y,ao} and the sextet {p^ai} among others. On the other hand 
under the U^(l) symmetry aU flavor multiplets are parity doubled, so the low-lying 
multiplets are the sextets {7r,ao} and {/?, ai} and the doublet {c,?^} among others. 

Ordinarily we would examine the multiplet degeneracy by measuring hadron 
masses. These are obtained by measuring hadron propagators between points widely 
separated in Euclidean time. At nonzero temperature, however, the Euclidean time 
dimension is limited, so we measure hadron propagation between points widely sep- 
arated in the spatial direction instead. Such propagators are often called "screening 
propagators" or "screening correlators". Eor example, from the generic meson prop- 
agator 

Gabir) = {HaiO)H,ir)), (48) 

where Ha{r) = t/j[r)Ta4^{r) with r = (x, j/, 2;, r), we may construct the static correlator 
by integrating over all but the z coordinate: 

Cab{z) = SabCa{z) = Sab j (It dxdyG a{x , y , z , t) (49) 

This correlator receives contributions from all possible states in a given channel. Thus 
the asymptotic behavior is given by 

oo 
n=0 

At large enough z the correlator falls off as 

Ca{z) ~ Paoe-""^"' (5f) 

where m^o is the lightest screening mass for the particular channel. These are the 
masses that may be used for exploring the multiplet structure. 

The screening propagators have a physical interpretation. Eor example, the pion 
screening state is the exchanged object that give rise to the Yukawa interaction be- 
tween two static nucleons in the thermal ensemble. 

The screening masses have been popular observables in nonzero temperature sim- 
ulations for some time."^^''*^''*"^''*'* An example showing the spectrum as a function of 
temperature is given in Eig. 9. The state labeled a in the flgure is really only the 
valence quark component, corresponding to the flrst term in Eq. (28). Let us call it 

39 

It is apparent that rapid changes take place in the screening spectrum at the 
phase transition. Can we tell from the multiplet structure which chiral symme- 
try is restored? To answer the question, we must take the chiral limit, i.e., take 
the quark masses to zero. The result appears to be consistent with the forma- 
tion of an SU(2) X SU(2) {7r,(T„} multiplet among others, suggesting restoration of 
SU(2) X SU(2). In these simulations the other {?y,ao} multiplet was not studied, so 
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Figure 9: Screening masses as a function of 6/(/^ (hence temperature).'*'* Results for 
four flavors of staggered quarks on a f6"^ X 8 lattice with ma = O.Of are shown in the 
left panel. The phase transition takes place at = 5.f5. The dot-dashed lines on 
the left mark from top to bottom, respectively, the zero temperature masses in units 
of Tc for the nucleon and /?, cr, and tt mesons. The dotted lines on the right indicate 
from top to bottom the screening masses expected from propagation of three and two 
free quarks. Results for a zero flavor (quenched) simulation on a 16"^ X 4 lattice are 
shown similarly in the right panel. The points at = 5.865 are from Ref. 43 and 

correspond to T = 1.5rc. The dot-dashed lines mark the zero temperature masses of 
the nucleon, /?, and tt. The free quark values differ in the two panels because of the 
different lattice sizes. 



we can't tell whether the Ua(1) symmetry is restored simultaneously. But a compli- 
cation arises. The valence part (Ty is not distinguishable from an ao in the continuum 
limit. '*° That would make the observed multiplet structure {7r,ao}, consistent with 
U^(l) restoration. However, it is argued in Ref. 39 that in the chiral limit a chiral 
selection rule eliminates the second (nonvalence) contribution to the a propagator, 
leaving us a multiplet structure consistent with the restoration of either symmetry. 
Further work is evidently needed to resolve the question. 

Shown on the high temperature side of the phase transition in Fig. 9 is the result 
expected from a continuum consisting of two and three free quarks. It is remarkable 
that most of the screening masses are in rough agreement with these values. In 
particular, the p and ai screening masses are close to the two-quark continuum and 



the nucleoli mass is close to the three- quark continuum. It is clear from Fig. 9 
that only the pion deviates significantly from the free-quark rule. Although such 
a result suggests deconfinement in the high temperature phase, the spectrum alone 
isn't decisive: charmonium has a mass close to twice the charm quark mass, but it is 
confined. 

The continuum value of the free quark screening mass is simply 



where the Lagrangian quark mass is niq. This result follows from the requirement 
that the quark field be antiperiodic in Euclidean time. Therefore the minimum time 
component of the quark and antiquark momentum is ttT and it propagates in the 
spatial direction with the stated effective mass. For sufficiently light quarks, the 
quark- antiquark threshold in the meson screening channel is 2TrT . For the baryon 
channel the three-quark threshold is SttT. These continuum values are modified on a 
finite lattice, as shown in the figure.'*'* 

Are these high temperature screening states merely a multiquark continuum? We 
return to this question in Sec. 9. 

6. Thermodynamics with Wilson Fermions 

To be confident that numerical simulations accurately represent continuum QCD 
it is essential that we demonstrate that our answers are independent of the fermion 
scheme. Unfortunately, thermodynamic simulations with Wilson fermions are not 
sufficiently developed at present to make a credible comparison with the staggered 
scheme. The fundamental difficulty is that we are dealing with a phase transition 
associated with the restoration of a spontaneously broken chiral symmetry, but the 
Wilson scheme builds in an explicit breaking of this symmetry, which goes away only 
in the continuum limit. 

For a given lattice dimension the staggered fermion parameter space consists of 
the gauge coupling and the quark masses ma. The chiral limit is reached at 

zero quark mass. In the Wilson fermion scheme each quark mass is replaced by a 
hopping parameter k. The chiral limit is not known a priori, but in simulations at 
low temperature and reasonably high values of 6/(/^, it is found that a number of 
indicators of chiral symmetry, e.g. a vanishing zero temperature pion mass and a 
vanishing current quark mass, hold at least with some consistency along a "chiral 
curve" in the parameter space 



among the indicators of chiral symmetry, so they must be checked in simulations. 
The thermal crossover occurs along another line depending on Nf 
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Figure 10: Phase diagram for two flavors of Wilson fermions from Ref. 20 as a function 
of the Wilson hopping parameter n and coupling fi = G/g^. The top curve k.^ locates 
the chiral line. The bottom curve gives the thermal crossover k.^ for = A and 
the middle curve gives the thermal crossover k.^ for = 6. In each case the low 
temperature phase is below the Kf curve and the high temperature phase is above. 
The arrows indicate the intersection of the thermal and chiral curves in each case as 
reported by the Tsukuba group. Data are from Refs. 49,46,47,48,45. 

The chiral line and thermal lines for A^t = 4 and Nf = 6 are indicated in Fig. 10. As 
the chiral line is approached from the low k, side, the quark mass is decreasing toward 
its physical value. We are interested in studying the crossover or phase transition in 
a region of small (physical) quark mass, so close to the chiral line. To avoid lattice 
artifacts, we would also like to work at a small lattice spacing. Indeed, evidence for 
lattice artifacts in two flavor simulations at Nf = 6 was recently reported by the MILC 
collaboration, which found a flrst order, possibly bulk, phase transition in simulations 
at /€ = 0.17, 0.18, and 0.19 and in the range 4.8 — 5.2. in close proximity to the 
thermal crossover.'*'' Thus one might infer that it would be preferable to work with 
values of larger than about 5. 

An unfortunate feature of the phase diagram, evident in Fig. 10, is that with 
A^f = 4, before we reach the chiral limit at large 6/(/^, we must flrst cross to the high 
temperature phase.'*® In other words, to study the crossover at small quark mass, we 
must work near the intersection of the thermal and chiral lines. The Tsukuba group 




discovered that for = A the intersection occurs at very small values of 6/(/^, usually 
considered to be in the strong coupling regime. At A^t = 6 the intersection occurs at 
slightly larger 6/(/^, as we can see in Fig. 10."*^ 
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Figure 11: From Ref. 50. Chiral condensate vs. n and various N-t for two flavors of 
Wilson fermions at 6/(/^ = 5.3. Dat from Refs. 45,47,50. 

At A^t = 8, not shown in the figure, another shift in the right direction is found. ^° 
A simulation at A^^ = 8 was done at 6/(/^ = 5.3 over a range of n up to Hc ~ 0.168, 
allowing a comparison with previous results at A^t = 4 and 6. The thermal crossover, 
now shifted to Ht ~ 0.167, shows no evidence for the presumed lattice artifact seen at 
slightly larger /€ at A^t = 6. As illustrated for {^V'V') in Fig. 11, as A^^ is increased, bulk 
quantities appear to follow an envelope established by the zero temperature theory, 
breaking away at the crossover. Moreover {^V'V') appears to be decreasing immediately 
prior to the crossover, suggesting progress toward a low quark mass.^° 

So, although the thermal line continues to shift in the right direction at A^t = 8, it 
is still not enough. Apparently it is necessary to work at quite high A^t, if we want a 
simulation that is as close to the continuum limit as the staggered fermion simulations 
whose predictions we are trying to confirm."*^ Thus with the original Wilson action, 
we are forced to chose between a strong coupling simulation and risk encountering 
lattice artifacts or carry out an expensive simulation. Further progress with the 
Wilson scheme is likely to require working with an improved action. The Tsukuba 



group has adopted one such improvement with encouraging prehminary results. 

7. Structural Changes at the Crossover 

Dramatic changes take place at the thermal crossover. Hadrons grow and merge, 
resulting in an extended mixture of quarks and gluons. Understanding the crossover 
is vitally important to the development of models of hadronization. In this section we 
examine two indicators of structural change in the QCD ensemble: the "constituent 
quark free energy" and the baryon density induced by the introduction of a point 
color triplet. We also consider the heavy-quark potential in a pure gluon ensemble. 
Finally we mention briefly recent efforts to explain the crossover in terms of topological 
structures, namely instantions and monopoles. 
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Figure 12: Constituent quark free energy as defined in text with c adjusted to 0.5. 
Data are from the same sources as Fig. 4. 



7.1. Constituent Quark Free Energy 

When a point color triplet (a spinless test quark) is placed in the thermal ensemble 
at low temperature, we expect that a gluon cloud surrounds it and a light antiquark 
binds to it to form a color singlet. This effect is induced in the ensemble by the 
Polyakov loop operator. The expectation value of the Polyakov loop variable measures 



the change in the free energy of the thermal ensemble due to the introduction of a 
point spinless test quark. This free energy difference, Eq. (30) 

/(r,m,) = -riog(ReP/3)c/, (55) 

is a function of the temperature T and light quark mass nig. At low temperature, 
the point source is screened by a light antiquark, forming a sort of heavy-light me- 
son. The free energy /(T, niq) then consists of the energy of the screening cloud and 
the lattice-regulated ultraviolet-divergent self-energy of the point source. The former 
contribution might be called loosely the "constituent quark free energy". At high tem- 
perature, screening is accomplished through thermal fluctuations in the color fields. 
The point source self-energy still diverges in the same way as at low temperature, 
however. Thus if we could isolate the contribution to /(T, m^) from the screening 
cloud, we would obtain the constituent quark free energy as a function of temperature, 
which may be of some interest in constructing models of the quark plasma. 

To compute the self energy of an isolated point source requires introducing an 
infrared cutoff, which can be temperature dependent. This requirement introduces 
an element of arbitrariness in the definition. In simplest terms, we must decide where 
to put the knife when we dissect the heavy meson. 

For an initial stab at this analysis, we observe that the regulated ultraviolet di- 
vergence is proportional to I/a. Thus we have 

f{T,mg)=fcq{T,mg) + c/a (56) 

for some constant c. Since Nf = 1 /aT^ we solve for the constituent quark free energy 

/cq(r, nig) = -r(log(ReP/3) + cN,). (57) 

This expression determines the scaling of the Polyakov loop measurement as Nf is 
varied. If the dimensionless constant c is chosen properly, and if continuum scaling 
holds, Polyakov loop data at increasing should yield a universal function /cq(r, rUg). 

The treated Polyakov loop data for a wide range of Nf is plotted in Fig. 12. The 
constant c is adjusted by eye to achieve the rough scaling agreement shown. For this 
purpose only values for the lightest available quark mass from each data set are used. 
Although the quark mass values nig/T are not the same from one to the next in this 
figure, they are small [nig/T < O.I) and would be expected to contribute little (of the 
order 10 MeV) to the free energy. Thus one would expect only a small inconsistency 
from this variation. It is amusing that at the crossover, the free energy drops by 
about the 300 MeV expected in a constituent quark model with deconfinement at 
high temperature. 

7.2. Induced Baryon Number 

Useful insights into the structure of the plasma can be obtained by measuring the 
quark number density in the vicinity of the point test charge. From our discussion of 




Figure 13: Quark number density induced by a fixed quark at the origin as a function 
of distance from the origin at three temperatures. Curves are fits to a single screening 
mass. The total induced quark number Q is also shown. 

the constituent quark self energy, we observed that at low temperature the screening 
cloud very likely contains an antiquark. Two light quarks may also combine with 
the point charge to form a baryon-like color singlet, but (at least we would expect) 
with lower probability, since the Boltzmann factor should suppress this contribution. 
Crudely, the suppression should go as exp(— mconst/T) where mconst is the mass of 
the extra constituent quark. At crossover, this mass is nearly twice the temperature, 
which would give an order of magnitude suppression. However, things are changing 
rapidly at the crossover. We learn more about the effect by measuring the dynamical 
quark number density Pq{r) in the presence of the test charge as a function of the 
distance from the test charge. The total induced quark number 

Q = j d\p,{r) (58) 

gives another measure of the screening cloud. 

Results of simulations with two light quark fiavors at three temperatures, O.TSTc, 
Tc, and l.bTc are shown in Fig. 13.^^. We see that as the temperature is increased, the 
correlation vanishes, just as one would expect if the plasma undergoes a crossover from 
a confined phase to a weakly correlated high temperature phase. The total induced 
charge, tabulated with the temperature in the figure, provides further information. 



Notice that at O.TSTc, the induced quark number is considerably different from — f . As 
we noted above, departures from — f are caused by diquark screening. A simple model 
based on the lowest S-wave baryons and mesons that could be formed in the screening 
process predicts Q = —0.81 at this temperature, considerably closer to —1.^^ To get 
the observed value would therefore require many more baryons than just the low lying 
S-wave states. Thus there appears to be an anomalously large baryonic component. 
If this effect is corroborated in further simulations, it might suggest a proliferation of 
baryons (and antibaryons) in the hadronic plasma as the crossover is approached from 
below. ^"^ Many years ago, Hagedorn proposed that an exponentially growing density 
of hadronic states would lead to a limiting temperature of about 160 it 10 MeV.^'* Of 
course, we now understand that there is no limiting temperature, since the entropy of 
hadronic states is limited by the entropy of the underlying quark and gluon degrees 
of freedom. However, the tendency for the density of hadronic states, particularly 
baryons, to grow rapidly with mass may lead to an entropy-induced proliferation of 
these states just below the crossover. 
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Figure 14: Heavy quark potential in an SU(3) gluon ensemble at various temperatures 
relative to the deconfinement temperature Tc . From Ref. 55. 



7.3. Heavy Quark Potential 

The correlation of two Polyakov loops (31) 

exp[-y(r,r)/r] = {P{0)P\r)/9)u 



(59) 



measures the change in free energy V{r^ T) of the thermal ensemble brought about by 
introducing a static quark and antiquark pair. It is a measure of the quark-antiquark 
potential (free-energy) as a function of separation distance r and temperature T . If 
dynamical quarks are omitted, at low temperature we expect the potential to grow 
linearly with distance, reflecting confinement. At high temperature deconfinement 
predicts a potential that becomes asymptotically constant at large r. This expectation 
is borne out in Fig. 14. With light quarks included the heavy quark is screened and 
we expect the potential to approach a constant at any temperature. 

7.^. Topological Structures 

The effort to develop phenomenological models of QCD, particularly to explain 
confinement, has frequently employed topological structures, particularly instantons 
and monopoles. These structures are to be found in lattice configurations. They 
may play an important role in the high temperature crossover, and they may reveal 
themselves as nonperturbative contributions to the equation of state. Although our 
understanding of these structures is still rudimentary, we mention them here for 
completeness. 

7 . 1 .Instantons 

Instantons are typically exposed in lattice configurations after a smoothing or 
cooling treatment to eliminate high frequency noise. A measure of the importance 
of instanton configurations is the topological susceptibility. Folklore says that at high 
temperature instantons should be suppressed, the axial gauge anomaly should disap- 
pear, and the U^(f ) chiral symmetry should be restored."^® Indeed, it was found some 
time ago in pure SU(N) gauge simulations that the susceptibility, measured after a 
cooling treatment, falls dramatically at the deconfinement transition.^'' DiGiacomo et 
al. cautioned, however, that the effect depends on way the susceptibility is defined.^® 
Thus if we want to correlate a falling topological susceptibility with the suppression 
of the gauge anomaly and the restoration of U^(I), it is important that we find a 
consistent, renormalizable definition. 

If instantons play a role in chiral symmetry breaking, then we should see a correla- 
tion between topological susceptibility and the chiral condensate. Recent progress in 
establishing this connection was reported by Teper.^^ In the Stony Brook instanton 
program, instantons take center stage in controlling the restoration of chiral sym- 
metry and in determining the screening spectrum. In recent work by Ilgenfritz and 
Shuryak and by Schafer, Shuryak, and Verbaarschot it is argued that in full QCD, as 
the chiral phase transition is approached, the light fermion determinant induces an at- 
tractive interaction between instantons and anti-instantons.^° The resulting molecules 
are predicted to predominate over solitary instantons and anti-instantons. The strong 



correlation drives the chiral phase transition. Schafer, Shuryak, and Verbaarschot 
compute hadronic screening masses in the model and find a spectrum in qualitative 
agreement with results from lattice simulations. It would be interesting to test their 
proposals further in lattice simulations. 

7.4.2.Monopoles 

A study of the thermal behavior of QCD may provide insight into the mechanism 
of confinement. An old 't Hooft — Mandelstam model characterizes the confining QCD 
vacuum as a dual superconductor, with an electric Meissner effect and a condensate 
of color magnetic monopoles.^^ Some years ago Schierholz et al and Kronfeld et al 
explored the association between confinement and the presence of monopole currents 
in Yang-Mills theories. Interest has revived recently.^"^'^'* 

To identify monopole currents in a nonabelian gauge theory it is necessary to 
carry out a U(I) projection of the gauge links. This is done by first fixing a suitable 
gauge. Popular choices include maximal Abelian gauge and a variety of "unitary" 
gauges, one of which involves diagonalizing the product of gauge links forming the 
Polyakov loop, inviting the conclusion that monopoles are responsible for confinement. 
A "monopole" contribution is then extracted from the resulting U(I) gauge field 
following the procedure of DeGrand and Toussaint.^^ 

In the new approach the string tension and other confinement features are com- 
puted using only the monopole contribution to the U(I) field. Good agreement is 
found with the full string tension computed in the conventional way. In the past year 
the Kanazawa group has also calculated the projected U(I)-monopole Polyakov-loop 
expectation value in both SU(2) and SU(3) Yang-Mills theory. The behavior of the 
monopole-projected Polyakov loop variable imitates the behavior of the conventional 
Polyakov loop variable. The similarity is seen in a variety of U(I) projection schemes. 

Although the results are promising, further work is needed, first, to demonstrate 
that the suitably defined abelian monopole currents survive the continuum limit^'' 
and, second, to find a suitable order parameter for monopole condensation.^® 

8. Bulk Properties of the Quark Gluon Plasma 

Three observables give important information about the characteristics of the 
quark plasma. The equation of state, giving energy density and pressure as a function 
of temperature, shows behavior expected in a rapid crossover from a phase dominated 
by hadrons to a phase dominated by quasi-free quarks and gluons. Departures from 
an ideal relativistic quark-gluon gas are particularly noticeable in the pressure for 
Tc < T < 2Tc. The baryon susceptibility, which measures the fiuctuation in baryon 



number, also shows behavior typical of deconfinement . 



8.1. Equation of State 

Phenomenological models of the high temperature phase require knowing its equa- 
tion of state and related bulk thermodynamic quantities, such as the rise in energy 
density at the crossover, the velocity of sound, and the contribution to the energy 
density from various species, including strange quarks. Particularly important is the 
determination of the peak in the specific heat near the phase transition for simulations 
with quarks. Simulations of the energy density and pressure of the quark plasma are 
costly because of a low signal-to-noise ratio. 

The earliest determinations of the energy density eiT^niq) and pressure piT^niq) 
as a function of temperature and quark mass"^ used the basic thermodynamic identities 



dF{V,T, 



m 



<T,mq)V = (60) 

p{T,mq) = — (61) 

In a lattice simulation each such derivative of the free energy involves a separate vari- 
ation of the spatial and temporal lattice constants a and aj, entailing a concomitant 
renormalization of the gauge coupling. Some years ago Karsch determined the per- 
turbative asymptotic variation of the gauge coupling with respect to the anisotropy 
parameter ^ = a^/a.^^ Unfortunately, present simulations are not in the perturbative 
asymptotic region. Although in principal the Karsch coefficients could be determined 
nonperturbatively from lattice simulations, this has not yet been done successfully.^^ 
Fortunately there is a different nonperturbative route to e and p. The "interaction 
measure" 

I = e-3p (62) 

is more easily determined, since it involves an isotropic variation of the lattice con- 
stant, requiring only the usual nonperturbative renormalization of the lattice coupling. 
The pressure, on the other hand, can be determined separately from the free energy 
density / = F/V = —p by integrating either of two relations™ 

, , dF(Q/g\m„,V) , , 

w„rt = '^''(;^;;;"^'' . (64) 

A vacuum subtraction is performed to give the pressure relative to the pressure of 
the nonperturbative vacuum: 
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Figure 15: From Ref. 25. (a) Fnergy density of the SU(3) gluon ensemble vs. tem- 
perature in units of the critical temperature T^, from simulations on 16"^ X 4, 32"^ X 6 
and 32"^ X 8 lattices. The dashed horizontal lines indicate the result expected for an 
ideal gas on these lattices. The solid horizontal line shows the continuum result. The 
vertical bar indicates the latent heat discontinuity, (b) Pressure in the SU(3) gluon 
ensemble vs. temperature. 

n faniq 

^ = iVM d(am') 

Of course the latter equation may be used only when dynamical quarks are present. 
Together with a determination of the interaction measure, this result can then be 
used to determine the energy density. 

These methods have been used in a nonperturbative determination of the energy 
density and pressure in SU(2) and SU(3) Yang-Mills theory^^'^^'''^ and in SU(3) with 
two flavors of dynamical staggered fermions.^® The pure glue theory can be simulated 
with higher precision than the theory with quarks included. It tells us in some 
detail how well QCD resembles a relativistic ideal gas in its bulk thermodynamic 
behavior. Figure 15 shows the result for the equation of state for the pure glue 
ensemble in SU(3). The curves are smooth having been derived by integrating a 
smooth interpolation of the plaquette measurements. 

The equation of state for a relativistic ideal gas for SU(N) color with Nj quark 
flavors is 

e = "^(^N'-I + InNj^^ (67) 
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(66) 



P = e/3 



(68) 



The continuum value must be corrected for the lattice discretization. This has been 
done in Fig. 15. It is clear that to a good approximation, for T > 2Tc^ as the contin- 
uum limit is approached, the equation of state resembles closely that of a relativistic 
ideal gas. However, the approach to the continuum Stefan-Boltzmann limit is slow. 
For T < 2Tc^ the pressure departs significantly from the free gas value and the en- 
ergy density differs by about 10-20%. One might have expected a strong deviation 
in the pressure, since mechanical stability requires that it be continuous across a 
phase transition. It has been suggested that the low momentum gluon modes, those 
that are most affected by strong interactions, may be responsible for the deviation 
from ideality. Thus precision studies of the equation of state may give indirect evi- 
dence for nonperturbative effects, such as contributions from instanton or monopole 
configurations in the thermal ensemble. 
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Figure 16: Fnergy density (upper curve) and three times the pressure (lower curve) 
vs. temperature (scale based on the rho meson mass) for two light quark fiavors 
{niq/T = 0.1) in the staggered fermion scheme from Ref. 26. The bursts give the 
pressure extrapolated to zero quark mass 

When quarks are included there is no evidence for a bona fide phase transition 
at nonzero quark mass. Nonetheless, there is a steep rise in the energy density at 
the temperature associated with the largest slope in the Polyakov loop and {^V'V'), as 
seen in Fig. 16. The transition region is remarkably sharp — of order 20 MeV. In a 
cooling quark plasma such a strong crossover could cause a momentary slowing in the 



expansion of the plasma as the quarks and gluons reorganized themselves into more 
compact hadrons. 
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Figure 17: Singlet baryon susceptibility vs. temperature for three lattice sizes and 
for two flavors of staggered fermions from Ref. 24. The horizontal line indicates the 
continuum ideal gas prediction. 



8.2. Baryon Susceptibility 

Most QCD simulations are done in the grand canonical ensemble (with respect to 
baryon number) at zero baryon chemical potential. As we have mentioned, technical 
difficulties have so far thwarted attempts at successful simulations at nonzero chemical 
potential, at least for credibly large lattice volumes. ^"^ This situation is unfortunate, 
since in heavy ion collisions one expects regions of nonzero baryon density. In the 
grand canonical ensemble at zero chemical potential, the average baryon density is 
zero. However, there are fluctuations. The baryon susceptibility quantifies these 
fiuctuations: 

dB{iJ.) 



Xs 



dp. 



(69) 



fM = 



At low temperature a fluctuation in baryon number requires the addition of a nu- 
cleon or antinucleon to the statistical ensemble. Thus we expect fluctuations to be 
suppressed by a Boltzmann factor exp(— mjv/r), where mjv is the nucleon mass. If 
the high temperature phase is characterized by an ideal gas of quarks and gluons, 
fluctuations are then controlled by the free energy of a quark or antiquark. The ideal 
gas (continuum) result is 

Xs = NfT' (70) 

The susceptibility can be measured separately for each quark flavor. For example, 
with two quark flavors, we may deflne chemical potentials //„ and fi^ for the up and 
down quarks, respectively. Then besides the singlet susceptibility we have deflned 
above, we have a nonsinglet susceptibility: 

_ dB dB 
"^"^ djiu djid 

which measures fluctuations in isospin. 

The baryon susceptibility has been measured in a few recent simulations.®^ Results 
are shown in Fig. 17. We see a rapid rise in baryon susceptibility at the crossover, 
as would be expected from deconflnement. There are signiflcant lattice discretization 
corrections for small Nf that gradually disappear on larger lattices. After allowing 
for these corrections, we see that this quantity agrees quite well with the ideal gas 
prediction in the high temperature phase. 



(71) 



9. Correlations and Confinement in the Quark Gluon Plasma 

We have already discussed static screening correlators in conjunction with chiral 
symmetry restoration. In this section we return to hadron correlators and examine 
the remarkable difference between their imaginary timelike ( "temporal" ) and spacelike 
behavior. Finally, we discuss the conflning features of the high temperature spacelike 
correlators. 



9.1. Hadron Propagation at Imaginary Time 

The Fourier transform of the screening propagator (48) is the temperature Green's 
function G'(cu,p,r): 

6abGa{iO,p,T) = j dT(Px{Ha{0)Hb{T,X.)) exp{-iijjT - ip ■ X.) (72) 

The static screening correlator (49) is then 

Ca{z) = ^ [ dp^exp{ip^z)Gaiid = 0,p^ = 0,py = 0,p^,T). (73) 



The screening masses locate the poles in the temperature Green's function in p^, at 
inian- Rather than integrating out the transverse coordinates, more generally, we 
could consider measuring the screening spectrum at nonzero cu, p^^ and py to get 
more information about screening. However, most studies set these variables to zero. 
Thus the static screening correlator is sensitive to poles in the temperature Green's 
function at zero frequency and imaginary wavenumber. 

An entirely different domain of the temperature Green's function can be reached 
by measuring the corresponding temporal correlator 

C*,(t) = J d'x{H^{0)H,{T,^))exp{-tp-^), (74) 

which is related to the temperature Green's function through 

CO 

Clbir) = TJ2 exp{tio.,T)Gab{iOn, P = 0, T) (75) 

n=0 

for Matsubara frequency = 2TrnT. This temporal correlation function is periodic 
in the imaginary time variable r with period 1/T. The r dependence is controlled 
by poles (and the continuum states) in io. The spectral decomposition of the Green's 
function is given by 

Gab{0Jn,P,T) = : 76 

J-oo ZTT lUn — OJ 

The spectral density /^^^(a;', p, T) is extremely interesting, since it is controls real 
time excitations of the plasma. Can it be measured on the lattice? In principle, it 
can, but in practice, it is extremely difficult. The problem is that in an imaginary 
time simulation, we measure the Green's function only at the discrete Matsubara 
frequencies cOn = 2TrnT. In fact, because our lattices have a finite A^t, we know the 
Green's function only for a finite set n = 0,. .. ,A^t — I. The mathematical prob- 
lem, then, is to carry out an analytic continuation from the finite discrete set to the 
whole complex plane. To make matters worse, the continuation is based on data 
with statistical uncertainties. To proceed, therefore, we must introduce assumptions 
about the form of the spectral density. In numerical simulations in condensed mat- 
ter physics, it is possible to progress by making a "maximum entropy" assumption 
about the form of the spectral density and by collecting data for Nf 50 or more.''^ 
Present thermodynamic simulations in lattice QCD are very far from approaching 
this standard. 

Despite the difficulties in extracting real-time spectral information from temporal 
correlators, the Bielefeld group has shown that they can be used to obtain interesting 
qualitative information about the quark plasma. ''"^ They compare temporal correlators 
with what would be expected if the correlations were dominated by free quark and 
antiquark propagation. (See Fig. 18.) The results show that at low temperature. 




Figure 18: From Ref. 73 The temporal correlators for the staggered fermion pseu- 
doscalar (filled circles) and vector channels (filled squares) at /3 = 5.1 (a), 5.3 (b) and 
6.5 (c), shown as a function of tT. In (c) the free quark/antiquark results for the 
pseudoscalar (solid line) and vector (dashed line) correlators are also shown. 

the correlator has a form that can be fit approximately to a single low-lying mesonic 
state, whereas at high temperature, the correlator has the form roughly expected of 
a free quark/antiquark pair. 

9.2. Screening Wave Functions 

We have seen that bulk quantities, such as the energy density, pressure, baryon 
susceptibility, constituent quark free energy, and temporal correlators all behave at 
sufficiently high temperature as though the quark-gluon plasma were a relativistic 
ideal gas. Lest we begin to believe too strongly in an ideal gas description of the 
plasma, we now revisit the screening states. Recall that the pion screening state is the 
exchanged object that give rise to the Yukawa interaction between two static nucleons 
in the plasma. It is predominantly a quark-antiquark state. At low temperature, 
confinement binds the quark and antiquark to produce the traditional Yukawa pion 
field. At high temperature, one might expect deconfinement would dissolve the pion. 
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Figure 19: The pion wave function at T = and the screening pion correlation 
function at T ^ 1.5 X Tc from Ref. 74. 

leaving only a qq continuum. To explore the structure of this state, one can measure 
its "wavefunction" . 

How is a wavefunction measured in a lattice simulation? At zero temperature the 
wavefunction for the lightest state in a given channel can be obtained by measuring 
the correlation function 

V-lx) oc J2{0{t = 0) V(y, T)rV'(y + x, r)) (77) 
y 

The operator O creates a quark and antiquark at t = and F is the appropriate 
Dirac matrix for the desired meson. The Fuclidean time separation r should be made 
large enough so that contributions from excited states die off and the result becomes 
independent of r apart from an overall normalization factor. This correlation function 
is not gauge invariant, so requires gauge fixing. Coulomb gauge is a popular choice. 
With F = 75 we obtain the pion screening state. 

The measurement of the wavefunction of the screening state is accomplished in 
precisely the same manner, after interchanging the roles of the imaginary time coor- 
dinate and one of the spatial coordinates-say z. Coulomb gauge is then defined with 
respect to a fixed value of z. 

A comparison of the wavefunction measured at zero temperature and high tem- 
perature is given in Fig. 19 J'^ The "zero" temperature lattices of size 16"^ X 24 were 



generated with two flavors of staggered quarks at 6/(/^ = 5.445 and niq = 0.025. The 
high temperature wavefunctions were generated on lattices of size 16^ X 24 X 4 at the 
same mass, couphng, and flavor number. The high temperature lattice parameters 
correspond to a temperature of approximately 1.5rc. 

The low and high temperature results are strikingly similar, suggesting strong 
correlations at high temperature. The source operator O in these simulations was 
an uncorrelated product of "wall" operators. Thus any correlation must arise from 
the interaction of the quark and antiquark. The ideal gas result for this observable 
would give an uncorrelated result independent of separation r. Results for other 
channels, such as the p meson and nucleon also show very little change from low to 
high temperature.'''*'''^ 

9.3. Dimensional Reduction and Confinement at High Temperature 

The mysterious results for the screening wavefunctions can be explained quite 
simply from an analysis of the Euclidean path integral.'''' In a Euclidean space all 
directions are placed on an equal footing. Thus we could just as well call one of the 
spatial directions, say the z direction, our Euclidean time axis r', and call the old 
periodic Euclidean time direction a new periodic direction z' . Erom this point of view 
the partition function for nonzero temperature QCD (call it QCDj for "thermal") is 
equivalent to the partition function for zero temperature QCD in a three-dimensional 
space with one of the three dimensions periodic (call it QCD^ for "compact"). Such 
a variant of QCD is conflning, as we shall see. The screening states of QCDj are 
just the conflned hadronic states of QCDc and the screening wavefunctions reflect 
this conflnement. As we have observed, the quarks in QCD^ have an effective mass 
that grows asymptotically with temperature as ttT . At very high temperatures, the 
periodic coordinate z' becomes so compact, we may speak of a dimensionally reduced 
theory.''® In this world all quark states behave like conflned heavy quark states. Thus 
at high temperature it is natural that screening masses are quantized in multiples of 
ttT according to their quark content. 

If QCDc is conflning, it should be possible to measure the spatial string tension 
in the theory. This has been done to quite high accuracy by the Bielefeld group. 
The result is as good a demonstration of conflnement in the compactifled theory 
as has been done for the zero temperature four-dimensional theory.''^ A recent high 
precision result is shown in Eig. 20. Eor temperatures T > 2Tc the spatial string 
tension deduced from these measurements is proportional to g^{T)T, as expected 
from dimensional reduction arguments. 

Inasmuch as high temperature meson and baryon screening states are heavy-quark 
bound states, they lend themselves to a nonrelativistic treatment. The Stony Brook 
group has explored models of such states.''^ 

As we have mentioned, the compactiflcation of the Euclidean time dimension in 
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Figure 20: Heavy quark pseudopotential showing confinement in the spatial direction 
in SU(3) with quarks omitted. Resuhs from Ref. 76 on 32"^ X lattices, correspond 
to three temperatures in the "deconfined" phase with Nr = 2 (circles), 3 (triangles) 
and 4 (squares) at /3 = 6.0. 

QCD results in an effective three-dimensional theory. The four- vector gauge potential 
A^(a;) is reinterpreted as an adjoint Higgs scalar Aq[x) and a three-vector potential 
A'^[x) for i = f,2, 3. This observation leads to an intriguing question, important 
for theoretical insight and model building, at least: Is the resulting theory found 
in a Higgs phase or a confining phase? The former option entails a spontaneous 
breaking of the gauge symmetry. In the Higgs phase the infrared behavior of the 
theory is most likely controlled by 't Hooft-Polyakov monopoles,®° leading to U(l) 
confinement. Thus either scenario would result in spatial confinement. In a recent 
study of SU(2) pure Yang- Mills theory, Karkkainen et al found strong evidence that 
the theory chooses the confining phase-that is, the field Ao{x) does not develop a 



nonzero vacuum expectation value. 



10. Conclusions 

Through numerical simulations in the staggered fermion scheme, considerable 
progress has been made in the past several years toward establishing a consistent 
qualitative picture of the high temperature behavior of QCD. Dynamical quarks make 
a clear difference in the behavior of the thermal ensemble. Whether there is a phase 
transition or only a crossover at physical values of the quark masses has not been 
rigorously established, but it is likely that any first order phase transition has a small 
latent heat. The crossover temperature is approximately in the range Tc = f40 — f 60 
MeV. The transition from the confining regime to the plasma regime takes place over 
a relatively narrow range of temperatures (approximately 20 MeV). Viewed over this 
range, the effective latent heat is large. The plasma is well characterized in bulk 
as an ideal gas of quarks and gluons for T ^ 2Tc^ but it retains characteristics of 
confinement, revealed in long-range screening phenomena, i.e. over distances greater 
than 0(1/ g^T). Strong deviations from ideality occur below about 2Tc. 

Some gaps remain in our qualitative understanding. The phase structure of QCD 
with the strange quark included ("2-|-f fiavors") has not been explored as thoroughly 
as the two fiavor theory. Particularly interesting is the study of critical behavior and 
soft modes in the multiparameter space of light quark masses. Thermodynamics in 
the Wilson fermion scheme has yet to establish itself as a contender. Since we need to 
be confident that the staggered fermion scheme is not misleading us, we must develop 
better methods for incorporating Wilson quarks. Finally, we need to develop good 
phenomenological models of the crossover. Topological models offer an intriguing 
direction. Further work is needed to establish their credibility. 
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